###############################################################################################################
#This is the full replication file for Rozenas, Schutte, and Zhukov. 2017. "The Political Legacy of Violence: #
#The Long-Term Impact of Stalin’s Repression in Ukraine". Journal of Politics (forthcoming).                  #
###############################################################################################################

#This file replicates the two empirical analyses for the main paper,
#as well as the additional tests conducted for the online appendix. 

########################################
#Loading libraries and replication data#
########################################

rm(list=ls())

#if required packages are missing, the lines below will install them:
list.of.packages <- c("rdrobust","spikes","rdd","AER","multiwayvcov","maptools","sp","spdep","GISTools","geosphere","rgdal","gtools","car","gdata", "nlme","xtable","lmtest","stargazer","ggplot2","latticeExtra","mgcv")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages)}
lapply(list.of.packages, require, character.only = TRUE)
rm(list.of.packages, new.packages)

mydir <- "~/Dropbox/SovietData/Replication/" # change path to user's local directory
load(paste0(mydir,'replication.RData')) 

################################################################################
############################Map the study region################################
################################################################################

#############################################
## Figure 1 of main text (map)
#############################################

sp_oblasts <- as.character(unique(rayons2$ra_oblast))

spOblast1 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[1]))
spOblast2 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[2]))
spOblast3 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[3]))
spOblast4 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[4]))
spOblast5 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[5]))
spOblast6 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[6]))
spOblast7 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[7]))
spOblast8 <- gUnionCascaded(subset(rayons2,ra_oblast==sp_oblasts[8]))

#Remove inner boundaries
spOblast1 <- gBuffer(spOblast1,byid=F,-0.1)
spOblast2 <- gBuffer(spOblast2,byid=F,-0.1)
spOblast3 <- gBuffer(spOblast3,byid=F,-0.1)
spOblast4 <- gBuffer(spOblast4,byid=F,-0.1)
spOblast5 <- gBuffer(spOblast5,byid=F,-0.1)
spOblast6 <- gBuffer(spOblast6,byid=F,-0.1)
spOblast7 <- gBuffer(spOblast7,byid=F,-0.1)
spOblast8 <- gBuffer(spOblast8,byid=F,-0.1)

#Generate layers for oblasts
oblast11.layer <- list("sp.lines",spOblast1,col="black",lty=1,lwd=6,first=FALSE)
oblast22.layer <- list("sp.lines",spOblast2,col="black",lty=1,lwd=6,first=FALSE)
oblast33.layer <- list("sp.lines",spOblast3,col="black",lty=1,lwd=6,first=FALSE)
oblast44.layer <- list("sp.lines",spOblast4,col="black",lty=1,lwd=6,first=FALSE)
oblast55.layer <- list("sp.lines",spOblast5,col="black",lty=1,lwd=6,first=FALSE)
oblast66.layer <- list("sp.lines",spOblast6,col="black",lty=1,lwd=6,first=FALSE)
oblast77.layer <- list("sp.lines",spOblast7,col="black",lty=1,lwd=6,first=FALSE)
oblast88.layer <- list("sp.lines",spOblast8,col="black",lty=1,lwd=6,first=FALSE)

oblast1.layer <- list("sp.lines",spOblast1,col="white",lty=1,lwd=2,first=FALSE)
oblast2.layer <- list("sp.lines",spOblast2,col="white",lty=1,lwd=2,first=FALSE)
oblast3.layer <- list("sp.lines",spOblast3,col="white",lty=1,lwd=2,first=FALSE)
oblast4.layer <- list("sp.lines",spOblast4,col="white",lty=1,lwd=2,first=FALSE)
oblast5.layer <- list("sp.lines",spOblast5,col="white",lty=1,lwd=2,first=FALSE)
oblast6.layer <- list("sp.lines",spOblast6,col="white",lty=1,lwd=2,first=FALSE)
oblast7.layer <- list("sp.lines",spOblast7,col="white",lty=1,lwd=2,first=FALSE)
oblast8.layer <- list("sp.lines",spOblast8,col="white",lty=1,lwd=2,first=FALSE)

#...and rayons
rayon.layer <- list("sp.lines",rayons2,col="black",lwd=1,first=FALSE)

#Finally, map it all 
spplot(rayons2,zcol="po_r14parl_mar",
       sp.layout=list(rayon.layer,
                      oblast11.layer,oblast22.layer,oblast33.layer,
                      oblast44.layer,oblast55.layer,oblast66.layer,
                      oblast77.layer,oblast88.layer,
                      oblast1.layer,oblast2.layer,oblast3.layer,
                      oblast4.layer,oblast5.layer,oblast6.layer,
                      oblast7.layer,oblast8.layer#,
                      #oblast11.layer,oblast22.layer,oblast33.layer,
                      #oblast44.layer,ooblast55.layer,oblast66.layer,
                      #oblast77.layer,oblast88.layer
       ),
       col.regions=gray.colors(20),main="'Pro-Russian' vote margin in the 2014 \n parliamentary elections")

spplot(rayons2,zcol="ev_deport",
       sp.layout=list(rayon.layer,
                      oblast11.layer,oblast22.layer,oblast33.layer,
                      oblast44.layer,oblast55.layer,oblast66.layer,
                      oblast77.layer,oblast88.layer,
                      oblast1.layer,oblast2.layer,oblast3.layer,
                      oblast4.layer,oblast5.layer,oblast6.layer,
                      oblast7.layer,oblast8.layer#,
       ),
       col.regions=gray.colors(20),main="Deportations during the Ukrainian \n Nationalist insurgency")


#############################################
## Figure S1 of Appendix (residualized map)
#############################################

#Moving to residualized maps
r_map1 <- lm(po_r14parl_mar ~ as.factor(ra_oblast),data = rayons2)
rayons2$res_po_r14parl_mar <- NA
rayons2$res_po_r14parl_mar[!is.na(rayons2$po_r14parl_mar)] <- as.numeric(r_map1$residuals)

spplot(rayons2,zcol="res_po_r14parl_mar",
       sp.layout=list(rayon.layer,oblast11.layer,oblast22.layer,oblast33.layer,
                      oblast44.layer,oblast55.layer,oblast66.layer,
                      oblast77.layer,oblast88.layer,
                      oblast1.layer,oblast2.layer,oblast3.layer,
                      oblast4.layer,oblast5.layer,oblast6.layer,
                      oblast7.layer,oblast8.layer#,
       ),
       col.regions=gray.colors(20),main="'Pro-Russian' vote margin in the 2014 \n parliamentary elections (residualized)")

r_map2 <- lm(ev_deport ~ as.factor(ra_oblast),data = rayons2)
rayons2$res_ev_deport <- NA
rayons2$res_ev_deport[!is.na(rayons2$ev_deport)] <- as.numeric(r_map2$residuals)

spplot(rayons2,zcol="res_ev_deport",
       sp.layout=list(rayon.layer,oblast11.layer,oblast22.layer,oblast33.layer,
                      oblast44.layer,oblast55.layer,oblast66.layer,
                      oblast77.layer,oblast88.layer,
                      oblast1.layer,oblast2.layer,oblast3.layer,
                      oblast4.layer,oblast5.layer,oblast6.layer,
                      oblast7.layer,oblast8.layer#,
       ),
       col.regions=gray.colors(20),main="Deportations during the Ukrainian \n Nationalist insurgency (residualized)")

################################################################################
#########################Instrumental Variable Analysis#########################
################################################################################

#Copy replication data and geographic information into different subsets
subz  <- merged_data
subz. <- rayons2

#############################################
## Table 1 of main text (summary stats)
#############################################

sumstat <- function(data,varz,labs=NULL){
  sum.mat <- data.frame(Obs=apply(data[,varz],2,function(x){sum(!is.na(x))}),
                        Mean=apply(data[,varz],2,function(x){mean(x,na.rm=T)}),
                        Median=apply(data[,varz],2,function(x){median(x,na.rm=T)}),
                        SD=apply(data[,varz],2,function(x){sd(x,na.rm=T)}),
                        Min=apply(data[,varz],2,function(x){min(x,na.rm=T)}),
                        Max=apply(data[,varz],2,function(x){max(x,na.rm=T)}))
  if(length(labs)==0){labs <- varz}
  row.names(sum.mat) <- labs
  return(sum.mat)
}

sum.tab <- round(sumstat(merged_data,varz=c("ev_deport","ra_ptsn_ctrl","ra_ptsn_ops","po_r14parl_mar","po_r14pres_mar","po_r12parl_mar","po_r10pres_mar","po_r07parl_mar","po_r06parl_mar","po_r04pres_mar","ra_near_railrd","ra_forest","ev_reb_any","ra_agro2","ra_urban","ra_occupation_duration","ra_rulang2_1931")),4)
var.labz <- c("Deported individuals","Soviet partisan control","Soviet partisan operations","`Pro-Russian' margin (parl. 2014)","`Pro-Russian' margin (pres. 2014)","`Pro-Russian' margin (parl. 2012)","`Pro-Russian' margin (pres. 2010)","`Pro-Russian' margin (parl. 2007)","`Pro-Russian' margin (parl. 2006)","`Pro-Russian' margin (pres. 2004)","Distance to railway (km)","Fraction forested","Intensity of rebel violence","Percent arable land","Urbanization","Days under German occupation","Russian-speaking in 1931") 
row.names(sum.tab) <- var.labz; rm(var.labz)
xtable(sum.tab)


#############################################
## Anecdotes
## KS test statistics reported in footnotes 5 and 6.
#############################################

# Deportations and railroads

X <- rayons2[subz.,]

mean(X$ev_deport[X$ra_near_railrd<mean(X$ra_near_railrd)])
mean(X$ev_deport[X$ra_near_railrd>mean(X$ra_near_railrd)])
ks.test(X$ev_deport[X$ra_near_railrd<mean(X$ra_near_railrd)],X$ev_deport[X$ra_near_railrd>mean(X$ra_near_railrd)])

mean(X$ev_reb_any[X$ra_near_railrd<mean(X$ra_near_railrd)])
mean(X$ev_reb_any[X$ra_near_railrd>mean(X$ra_near_railrd)])
ks.test(X$ev_reb_any[X$ra_near_railrd<mean(X$ra_near_railrd)],X$ev_reb_any[X$ra_near_railrd>mean(X$ra_near_railrd)])

mean(X$ev_gov_any[X$ra_near_railrd<mean(X$ra_near_railrd)])
mean(X$ev_gov_any[X$ra_near_railrd>mean(X$ra_near_railrd)])
ks.test(X$ev_gov_any[X$ra_near_railrd<mean(X$ra_near_railrd)],X$ev_gov_any[X$ra_near_railrd>mean(X$ra_near_railrd)])


#######################################################
#######################################################
## MAIN IV RESULTS
## These models are shown in table 2 of the main text
## and tables 1 to 5 in the supplementary information
#######################################################
#######################################################

names(rayons2)

## Regression tables
cv <- "ev_reb_any + ra_agro2 + ra_occupation_duration + ra_urban"
pro.r <- c("po_r14parl_mar","po_r14pres_mar","po_r12parl_mar","po_r10pres_mar","po_r07parl_mar","po_r06parl_mar","po_r04pres_mar")

y.vec <- pro.r
x.vec <- c("ev_deport","ra_ptsn_ops")
z.vec <- c("ra_near_railrd","ra_forest")
X <- rayons2[subz.,]
cv.vec <- strsplit(cv," \\+ ")[[1]]
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
yx.vec <- data.frame(y.vec=rep(y.vec,length(x.vec)),x.vec=rep(x.vec,each=length(y.vec)))
for(j in 1:ncol(yx.vec)){yx.vec[,j] <- as.character(yx.vec[,j])}
yx.vec


# Model 1
i <- 1
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod1 <- ivreg(formz,data=X);  summary(mod1,diagnostics=T)
# first stage
mod1b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod1b,diagnostics=T)


# Model 2
i <- 2
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod2 <- ivreg(formz,data=X);  summary(mod2,diagnostics=T)
# first stage
mod2b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod2b,diagnostics=T)



# Model 3
i <- 3
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod3 <- ivreg(formz,data=X);  
summary(mod3,diagnostics=T)
# first stage
mod3b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod3b,diagnostics=T)

# Model 4
i <- 4
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod4 <- ivreg(formz,data=X);  
summary(mod4,diagnostics=T)
# first stage
mod4b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod4b,diagnostics=T)

# Model 5
i <- 5
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod5 <- ivreg(formz,data=X);  
summary(mod5,diagnostics=T)
# first stage
mod5b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod5b,diagnostics=T)

# Model 6
i <- 6
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod6 <- ivreg(formz,data=X);  
summary(mod6,diagnostics=T)
# first stage
mod6b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod6b,diagnostics=T)

# Model 7
i <- 7
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod7 <- ivreg(formz,data=X);  
summary(mod7,diagnostics=T)
# first stage
mod7b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod7b,diagnostics=T)

# Model 8
i <- 8
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod8 <- ivreg(formz,data=X);  
summary(mod8,diagnostics=T)
# first stage
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~ra_oblast + ra_ptsn_ctrl +ev_deport + ev_reb_any") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
ray.ME1 <- SpatialFiltering(paste("ra_ptsn_ctrl ~ ra_oblast + ",paste(z.vec,collapse="+")," + fitted(ray.ME2)") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
mod8b <- glm(paste("ra_ptsn_ctrl","~ as.factor(ra_oblast) +",paste(z.vec,collapse="+"),"+fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod8b,diagnostics=T)

# Model 9
i <- 9
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod9 <- ivreg(formz,data=X);  
summary(mod9,diagnostics=T)
# first stage
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~ra_oblast + ra_ptsn_ctrl +ev_deport + ev_reb_any") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
ray.ME1 <- SpatialFiltering(paste("ra_ptsn_ctrl ~ ra_oblast + ",paste(z.vec,collapse="+")," + fitted(ray.ME2)") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
mod9b <- glm(paste("ra_ptsn_ctrl","~ as.factor(ra_oblast) +",paste(z.vec,collapse="+"),"+fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod9b,diagnostics=T)

# Model 10
i <- 10
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod10 <- ivreg(formz,data=X);  
summary(mod10,diagnostics=T)
# first stage
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~ra_oblast + ra_ptsn_ctrl +ev_deport + ev_reb_any") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
ray.ME1 <- SpatialFiltering(paste("ra_ptsn_ctrl ~ ra_oblast + ",paste(z.vec,collapse="+")," + fitted(ray.ME2)") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
mod10b <- glm(paste("ra_ptsn_ctrl","~ as.factor(ra_oblast) +",paste(z.vec,collapse="+"),"+fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod10b,diagnostics=T)

# Model 11
i <- 11
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod11 <- ivreg(formz,data=X);  
summary(mod11,diagnostics=T)
# first stage
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~ra_oblast + ra_ptsn_ctrl +ev_deport + ev_reb_any") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
ray.ME1 <- SpatialFiltering(paste("ra_ptsn_ctrl ~ ra_oblast + ",paste(z.vec,collapse="+")," + fitted(ray.ME2)") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
mod11b <- glm(paste("ra_ptsn_ctrl","~ as.factor(ra_oblast) +",paste(z.vec,collapse="+"),"+fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod11b,diagnostics=T)

# Model 12
i <- 12
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod12 <- ivreg(formz,data=X);  
summary(mod12,diagnostics=T)
# first stage
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~ra_oblast + ra_ptsn_ctrl +ev_deport + ev_reb_any") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
ray.ME1 <- SpatialFiltering(paste("ra_ptsn_ctrl ~ ra_oblast + ",paste(z.vec,collapse="+")," + fitted(ray.ME2)") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
mod12b <- glm(paste("ra_ptsn_ctrl","~ as.factor(ra_oblast) +",paste(z.vec,collapse="+"),"+fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod12b,diagnostics=T)

# Model 13
i <- 13
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod13 <- ivreg(formz,data=X);  
summary(mod13,diagnostics=T)
# first stage
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~ra_oblast + ra_ptsn_ctrl +ev_deport + ev_reb_any") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
ray.ME1 <- SpatialFiltering(paste("ra_ptsn_ctrl ~ ra_oblast + ",paste(z.vec,collapse="+")," + fitted(ray.ME2)") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
mod13b <- glm(paste("ra_ptsn_ctrl","~ as.factor(ra_oblast) +",paste(z.vec,collapse="+"),"+fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod13b,diagnostics=T)

# Model 14
i <- 14
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod14 <- ivreg(formz,data=X);  
summary(mod14,diagnostics=T)
# first stage
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~ra_oblast + ra_ptsn_ctrl +ev_deport + ev_reb_any") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
ray.ME1 <- SpatialFiltering(paste("ra_ptsn_ctrl ~ ra_oblast + ",paste(z.vec,collapse="+")," + fitted(ray.ME2)") , data = X, nb = W_nb, style = "W",  ExactEV=T, zero.policy=T)
mod14b <- glm(paste("ra_ptsn_ctrl","~ as.factor(ra_oblast) +",paste(z.vec,collapse="+"),"+fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod14b,diagnostics=T)

#############################################
## Table 2 in the main text
## (also Table S1 in Appendix)
#############################################

# Make table (deport)
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
# 1st stage
stargazer(mod1b,mod2b,mod3b,mod4b,mod5b,mod6b,mod7b)
# Diagnostics
digz <- c("wi.test","wi.test.p","wu.test","wu.test.p","sa.test","sa.test.p","moran.st","moran.p")
digz.mat <- as.data.frame(matrix(NA,nrow=length(digz),ncol=7))
row.names(digz.mat) <- digz
colnames(digz.mat) <- 1:7
modz <- list(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
wi.test <- c(); wu.test <- c(); sa.test <- c(); wi.test.p <- c(); wu.test.p <- c(); sa.test.p <- c(); moran.st <- c(); moran.p <- c();
for(i in 1:7){
  X <- rayons2[subz.,]
  X <- X[!is.na(as.data.frame(X)[,strsplit(as.character(modz[[i]]$formula),"~")[[2]]]),]
  W_nb <- poly2nb(X, queen = T)
  W_listw <- nb2listw(W_nb, style = "W", zero.policy=T)
  digz.mat["wi.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,3]
  digz.mat["wu.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,3]
  digz.mat["sa.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,3]
  digz.mat["wi.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,4]
  digz.mat["wu.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,4]
  digz.mat["sa.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,4]
  digz.mat["moran.st",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$statistic
  digz.mat["moran.p",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$p.value
}
digz.st <- digz.mat[!grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- digz.mat[grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- ifelse(digz.p<.01,"**",ifelse(digz.p<.05,"*",ifelse(digz.p<.1,"'","")))
for(j in 1:ncol(digz.st)){digz.st[,j] <- as.character(round(digz.st[,j],3))}
for(i in 1:nrow(digz.st)){
  for(j in 1:ncol(digz.st)){
    digz.st[i,j] <- paste0(digz.st[i,j],digz.p[i,j])
  }}
xtable(digz.st)

#############################################
# Table S3 in Appendix (partisans)
#############################################

# Make table (2nd stage)
stargazer(mod8,mod9,mod10,mod11,mod12,mod13,mod14)
# 1st stage
stargazer(mod8b,mod9b,mod10b,mod11b,mod12b,mod13b,mod14b)
# Diagnostics
digz <- c("wi.test","wi.test.p","wu.test","wu.test.p","sa.test","sa.test.p","moran.st","moran.p")
digz.mat <- as.data.frame(matrix(NA,nrow=length(digz),ncol=7))
row.names(digz.mat) <- digz
colnames(digz.mat) <- 1:7
modz <- list(mod8,mod9,mod10,mod11,mod12,mod13,mod14)
wi.test <- c(); wu.test <- c(); sa.test <- c(); wi.test.p <- c(); wu.test.p <- c(); sa.test.p <- c(); moran.st <- c(); moran.p <- c();
for(i in 1:7){
  X <- rayons2[subz.,]
  X <- X[!is.na(as.data.frame(X)[,strsplit(as.character(modz[[i]]$formula),"~")[[2]]]),]
  W_nb <- poly2nb(X, queen = T)
  W_listw <- nb2listw(W_nb, style = "W", zero.policy=T)
  digz.mat["wi.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,3]
  digz.mat["wu.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,3]
  digz.mat["sa.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,3]
  digz.mat["wi.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,4]
  digz.mat["wu.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,4]
  digz.mat["sa.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,4]
  digz.mat["moran.st",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$statistic
  digz.mat["moran.p",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$p.value
}
digz.st <- digz.mat[!grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- digz.mat[grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- ifelse(digz.p<.01,"**",ifelse(digz.p<.05,"*",ifelse(digz.p<.1,"'","")))
for(j in 1:ncol(digz.st)){digz.st[,j] <- as.character(round(digz.st[,j],3))}
for(i in 1:nrow(digz.st)){
  for(j in 1:ncol(digz.st)){
    digz.st[i,j] <- paste0(digz.st[i,j],digz.p[i,j])
  }}
xtable(digz.st)

## Quantities of interest
X <- rayons2[subz.,]
sd(X$ra_near_railrd)
sd(X$ev_deport)*coefficients(mod1b)["ra_near_railrd"]
sd(X$po_r14parl_mar,na.rm=T)*coefficients(mod1)["ev_deport"]
sd(X$po_r12parl_mar,na.rm=T)*coefficients(mod3)["ev_deport"]
mean(X$ra_ptsn_ops)


#############################################
## Figure 2 in main text (coefficient plot)
#############################################

# Deportations

mod.list <- list(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
labs1 <- c(2014,2014,2012,2010,2007,2006,2004)
labs2 <- c("Parliament","President","Parliament","President","Parliament","Parliament","President")
xlims <- c(0.25,7.75)
ylims <- c(-.6,.6)
par(mar=c(2,2,.5,.5))
plot(0,0,xlim=xlims,ylim=ylims,xaxt="n",yaxt="n",xlab="",ylab="",bty="n",col=NA)
rect(xleft=xlims[1],xright=xlims[2],ybottom=ylims[1],ytop=ylims[2],col="grey90",border = NA)
segments(x0=xlims[1],x1=xlims[2],y0=pretty(ylims),y1=pretty(ylims),col="white",lwd=2)
segments(x0=xlims[1],x1=xlims[2],y0=pretty(ylims)-diff(pretty(ylims))/2,y1=pretty(ylims)-diff(pretty(ylims))/2,col="white",lwd=.5)
segments(x0=xlims[1],x1=xlims[2],y0=0,y1=0,col="gray40",lwd=3)
for(i in 1:length(mod.list)){
  submod <- mod.list[[i]]
  coefficients(submod)["ev_deport"]
  beta.hat <- summary(submod)$coefficients["ev_deport","Estimate"]
  sd.hat <- summary(submod)$coefficients["ev_deport","Std. Error"]
  rect(xleft=i-.25,xright=i+.25,ytop=beta.hat+1.96*sd.hat,ybottom=beta.hat-1.96*sd.hat,border=NA,col=rgb(.8,.2,.2,.5))
  segments(x0=i-.25,x1=i+.25,y0=beta.hat,y1=beta.hat)
}
axis(2,at=pretty(ylims),cex.axis=.8)
axis(1,at=1:7,labels = labs1,tick=FALSE,padj = -1.6,cex.axis=1.2)
axis(1,at=1:7,labels = labs2,tick=FALSE,padj = -.6,col.axis="grey",cex.axis=1)


#############################################
## Figure 4 in main text (coefficient plot)
############################################## 

# Partisans
mod.list <- list(mod8,mod9,mod10,mod11,mod12,mod13,mod14)
labs1 <- c(2014,2014,2012,2010,2007,2006,2004)
labs2 <- c("Parliament","President","Parliament","President","Parliament","Parliament","President")
xlims <- c(0.25,7.75)
ylims <- c(-.25,.25)
par(mar=c(2,2,.5,.5))
plot(0,0,xlim=xlims,ylim=ylims,xaxt="n",yaxt="n",xlab="",ylab="",bty="n",col=NA)
rect(xleft=xlims[1],xright=xlims[2],ybottom=ylims[1],ytop=ylims[2],col="grey90",border = NA)
segments(x0=xlims[1],x1=xlims[2],y0=pretty(ylims),y1=pretty(ylims),col="white",lwd=2)
segments(x0=xlims[1],x1=xlims[2],y0=pretty(ylims)-diff(pretty(ylims))/2,y1=pretty(ylims)-diff(pretty(ylims))/2,col="white",lwd=.5)
segments(x0=xlims[1],x1=xlims[2],y0=0,y1=0,col="gray40",lwd=3)
for(i in 1:length(mod.list)){
  submod <- mod.list[[i]]
  coefficients(submod)["ra_ptsn_ops"]
  beta.hat <- summary(submod)$coefficients["ra_ptsn_ops","Estimate"]
  sd.hat <- summary(submod)$coefficients["ra_ptsn_ops","Std. Error"]
  rect(xleft=i-.25,xright=i+.25,ytop=beta.hat+1.96*sd.hat,ybottom=beta.hat-1.96*sd.hat,border=NA,col=rgb(.2,.2,.8,.5))
  segments(x0=i-.25,x1=i+.25,y0=beta.hat,y1=beta.hat)
}
axis(2,at=pretty(ylims),cex.axis=.8)
axis(1,at=1:7,labels = labs1,tick=FALSE,padj = -1.6,cex.axis=1.2)
axis(1,at=1:7,labels = labs2,tick=FALSE,padj = -.6,col.axis="grey",cex.axis=1)



###########################################################
## Table S2 in Appendix (w/ 1931 census data as covariate)
###########################################################

cv <- "ev_reb_any + ra_agro2 + ra_occupation_duration + ra_urban + ra_rulang2_1931"
pro.r <- c("po_r14parl_mar","po_r14pres_mar","po_r12parl_mar","po_r10pres_mar","po_r07parl_mar","po_r06parl_mar","po_r04pres_mar")

names(X)

y.vec <- pro.r
x.vec <- c("ev_deport","ra_ptsn_ops")
z.vec <- c("ra_near_railrd","ra_forest")
X <- rayons2[subz.,]
cv.vec <- strsplit(cv," \\+ ")[[1]]
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
yx.vec <- data.frame(y.vec=rep(y.vec,length(x.vec)),x.vec=rep(x.vec,each=length(y.vec)))
for(j in 1:ncol(yx.vec)){yx.vec[,j] <- as.character(yx.vec[,j])}
yx.vec
summary(X[,cv.vec])

# Model 1
i <- 1
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod1 <- ivreg(formz,data=X);  
summary(mod1,diagnostics=T)

# Model 2
i <- 2
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod2 <- ivreg(formz,data=X);  
summary(mod2,diagnostics=T)


# Model 3
i <- 3
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod3 <- ivreg(formz,data=X);  
summary(mod3,diagnostics=T)

# Model 4
i <- 4
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod4 <- ivreg(formz,data=X);  
summary(mod4,diagnostics=T)

# Model 5
i <- 5
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod5 <- ivreg(formz,data=X);  
summary(mod5,diagnostics=T)

# Model 6
i <- 6
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod6 <- ivreg(formz,data=X);  
summary(mod6,diagnostics=T)

# Model 7
i <- 7
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod7 <- ivreg(formz,data=X);  
summary(mod7,diagnostics=T)


# Make table (deport)
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7)

# Diagnostics table (deport)
digz <- c("wi.test","wi.test.p","wu.test","wu.test.p","sa.test","sa.test.p","moran.st","moran.p")
digz.mat <- as.data.frame(matrix(NA,nrow=length(digz),ncol=7))
row.names(digz.mat) <- digz
colnames(digz.mat) <- 1:7
modz <- list(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
wi.test <- c(); wu.test <- c(); sa.test <- c(); wi.test.p <- c(); wu.test.p <- c(); sa.test.p <- c(); moran.st <- c(); moran.p <- c();
for(i in 1:7){
  X <- rayons2[subz.,]
  X <- X[!is.na(as.data.frame(X)[,strsplit(as.character(modz[[i]]$formula),"~")[[2]]]),]
  W_nb <- poly2nb(X, queen = T)
  W_listw <- nb2listw(W_nb, style = "W", zero.policy=T)
  digz.mat["wi.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,3]
  digz.mat["wu.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,3]
  digz.mat["sa.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,3]
  digz.mat["wi.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,4]
  digz.mat["wu.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,4]
  digz.mat["sa.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,4]
  #digz.mat["moran.st",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$statistic
  #digz.mat["moran.p",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$p.value
}
digz.st <- digz.mat[!grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- digz.mat[grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- ifelse(digz.p<.01,"**",ifelse(digz.p<.05,"*",ifelse(digz.p<.1,"'","")))
for(j in 1:ncol(digz.st)){digz.st[,j] <- as.character(round(digz.st[,j],3))}
for(i in 1:nrow(digz.st)){
  for(j in 1:ncol(digz.st)){
    digz.st[i,j] <- paste0(digz.st[i,j],digz.p[i,j])
  }}

xtable(digz.st)


#############################################################
## Table S4 of Appendix (effect of UPA violence on voting)
#############################################################

cv <- "ra_agro2 + ra_occupation_duration + ra_urban"
pro.r <- c("po_r14parl_mar","po_r14pres_mar","po_r12parl_mar","po_r10pres_mar","po_r07parl_mar","po_r06parl_mar","po_r04pres_mar")
y.vec <- pro.r
x.vec <- c("ev_reb_any")
z.vec <- c("ra_dist_to_ra","ra_forest")
X <- rayons2[subz.,]
cv.vec <- strsplit(cv," \\+ ")[[1]]
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
yx.vec <- data.frame(y.vec=rep(y.vec,each=length(x.vec)),x.vec=rep(x.vec,length(y.vec)))
for(j in 1:ncol(yx.vec)){yx.vec[,j] <- as.character(yx.vec[,j])}
yx.vec
summary(X[,cv.vec])

# Model 1
i <- 1
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod1 <- ivreg(formz,data=X);  
summary(mod1,diagnostics=T)
mod1b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod1b,diagnostics=T)

# Model 2
i <- 2
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod2 <- ivreg(formz,data=X);  
summary(mod2,diagnostics=T)
mod2b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod2b,diagnostics=T)


# Model 3
i <- 3
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod3 <- ivreg(formz,data=X);  
summary(mod3,diagnostics=T)
mod3b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod3b,diagnostics=T)

# Model 4
i <- 4
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod4 <- ivreg(formz,data=X);  
summary(mod4,diagnostics=T)
mod4b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod4b,diagnostics=T)

# Model 5
i <- 5
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod5 <- ivreg(formz,data=X);  
summary(mod5,diagnostics=T)
mod5b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod5b,diagnostics=T)

# Model 6
i <- 6
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod6 <- ivreg(formz,data=X);  
summary(mod6,diagnostics=T)
mod6b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod6b,diagnostics=T)

# Model 7
i <- 7
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod7 <- ivreg(formz,data=X);  
summary(mod7,diagnostics=T)
mod7b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod7b,diagnostics=T)


# Make table (deport)
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
# 1st stage
stargazer(mod1b,mod2b,mod3b,mod4b,mod5b,mod6b,mod7b)
# Diagnostics table (deport)
digz <- c("wi.test","wi.test.p","wu.test","wu.test.p","sa.test","sa.test.p","moran.st","moran.p")
digz.mat <- as.data.frame(matrix(NA,nrow=length(digz),ncol=7))
row.names(digz.mat) <- digz
colnames(digz.mat) <- 1:7
modz <- list(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
wi.test <- c(); wu.test <- c(); sa.test <- c(); wi.test.p <- c(); wu.test.p <- c(); sa.test.p <- c(); moran.st <- c(); moran.p <- c();
for(i in 1:7){
  X <- rayons2[subz.,]
  X <- X[!is.na(as.data.frame(X)[,strsplit(as.character(modz[[i]]$formula),"~")[[2]]]),]
  W_nb <- poly2nb(X, queen = T)
  W_listw <- nb2listw(W_nb, style = "W", zero.policy=T)
  digz.mat["wi.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,3]
  digz.mat["wu.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,3]
  digz.mat["sa.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,3]
  digz.mat["wi.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,4]
  digz.mat["wu.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,4]
  digz.mat["sa.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,4]
  digz.mat["moran.st",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$statistic
  digz.mat["moran.p",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$p.value
}
digz.st <- digz.mat[!grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- digz.mat[grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- ifelse(digz.p<.01,"**",ifelse(digz.p<.05,"*",ifelse(digz.p<.1,"'","")))
for(j in 1:ncol(digz.st)){digz.st[,j] <- as.character(round(digz.st[,j],3))}
for(i in 1:nrow(digz.st)){
  for(j in 1:ncol(digz.st)){
    digz.st[i,j] <- paste0(digz.st[i,j],digz.p[i,j])
  }}

xtable(digz.st)



#############################################################
## Table S5 of Appendix (effect of other Soviet violence)
#############################################################

cv <- "ra_agro2 + ra_occupation_duration + ra_urban"
pro.r <- c("po_r14parl_mar","po_r14pres_mar","po_r12parl_mar","po_r10pres_mar","po_r07parl_mar","po_r06parl_mar","po_r04pres_mar")
y.vec <- pro.r
x.vec <- c("ev_gov_any")
z.vec <- c("ra_dist_to_ra","ra_forest")
X <- rayons2[subz.,]
cv.vec <- strsplit(cv," \\+ ")[[1]]
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
yx.vec <- data.frame(y.vec=rep(y.vec,each=length(x.vec)),x.vec=rep(x.vec,length(y.vec)))
for(j in 1:ncol(yx.vec)){yx.vec[,j] <- as.character(yx.vec[,j])}
yx.vec
summary(X[,cv.vec])

# Model 1
i <- 1
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod1 <- ivreg(formz,data=X);  
summary(mod1,diagnostics=T)
mod1b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod1b,diagnostics=T)

# Model 2
i <- 2
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod2 <- ivreg(formz,data=X);  
summary(mod2,diagnostics=T)
mod2b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod2b,diagnostics=T)


# Model 3
i <- 3
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod3 <- ivreg(formz,data=X);  
summary(mod3,diagnostics=T)
mod3b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod3b,diagnostics=T)

# Model 4
i <- 4
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod4 <- ivreg(formz,data=X);  
summary(mod4,diagnostics=T)
mod4b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod4b,diagnostics=T)

# Model 5
i <- 5
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod5 <- ivreg(formz,data=X);  
summary(mod5,diagnostics=T)
mod5b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod5b,diagnostics=T)

# Model 6
i <- 6
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod6 <- ivreg(formz,data=X);  
summary(mod6,diagnostics=T)
mod6b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod6b,diagnostics=T)

# Model 7
i <- 7
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod7 <- ivreg(formz,data=X);  
summary(mod7,diagnostics=T)
mod7b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod7b,diagnostics=T)


# Make table
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
# 1st stage
stargazer(mod1b,mod2b,mod3b,mod4b,mod5b,mod6b,mod7b)
# Diagnostics table (deport)
digz <- c("wi.test","wi.test.p","wu.test","wu.test.p","sa.test","sa.test.p","moran.st","moran.p")
digz.mat <- as.data.frame(matrix(NA,nrow=length(digz),ncol=7))
row.names(digz.mat) <- digz
colnames(digz.mat) <- 1:7
modz <- list(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
wi.test <- c(); wu.test <- c(); sa.test <- c(); wi.test.p <- c(); wu.test.p <- c(); sa.test.p <- c(); moran.st <- c(); moran.p <- c();
for(i in 1:7){
  X <- rayons2[subz.,]
  X <- X[!is.na(as.data.frame(X)[,strsplit(as.character(modz[[i]]$formula),"~")[[2]]]),]
  W_nb <- poly2nb(X, queen = T)
  W_listw <- nb2listw(W_nb, style = "W", zero.policy=T)
  digz.mat["wi.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,3]
  digz.mat["wu.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,3]
  digz.mat["sa.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,3]
  digz.mat["wi.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,4]
  digz.mat["wu.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,4]
  digz.mat["sa.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,4]
  digz.mat["moran.st",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$statistic
  digz.mat["moran.p",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$p.value
}
digz.st <- digz.mat[!grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- digz.mat[grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- ifelse(digz.p<.01,"**",ifelse(digz.p<.05,"*",ifelse(digz.p<.1,"'","")))
for(j in 1:ncol(digz.st)){digz.st[,j] <- as.character(round(digz.st[,j],3))}
for(i in 1:nrow(digz.st)){
  for(j in 1:ncol(digz.st)){
    digz.st[i,j] <- paste0(digz.st[i,j],digz.p[i,j])
  }}

xtable(digz.st)



###########################################################
## Table S9 in Appendix (predictors of 2001 ethno-linguistic composition)
###########################################################

cv <- "ev_reb_any + ra_agro2 + ra_occupation_duration + ra_urban"
y.vec <- c("ra_ruminukr")
x.vec <- c("ev_deport","ra_ptsn_ops")
z.vec <- c("ra_near_railrd","ra_forest")
X <- rayons2[subz.,]
cv.vec <- strsplit(cv," \\+ ")[[1]]
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
yx.vec <- data.frame(y.vec=rep(y.vec,each=length(x.vec)),x.vec=rep(x.vec,length(y.vec)))
for(j in 1:ncol(yx.vec)){yx.vec[,j] <- as.character(yx.vec[,j])}
yx.vec
summary(X[,cv.vec])

names(X)

# Model 1
i <- 1
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod1 <- ivreg(formz,data=X);  
summary(mod1,diagnostics=T)

# Model 2
i <- 2
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod2 <- ivreg(formz,data=X);  
summary(mod2,diagnostics=T)


# Make table (deport)
stargazer(mod1,mod2)
# Diagnostics table (deport)
digz <- c("wi.test","wi.test.p","wu.test","wu.test.p","sa.test","sa.test.p","moran.st","moran.p")
digz.mat <- as.data.frame(matrix(NA,nrow=length(digz),ncol=2))
row.names(digz.mat) <- digz
colnames(digz.mat) <- 1:2
modz <- list(mod1,mod2)
wi.test <- c(); wu.test <- c(); sa.test <- c(); wi.test.p <- c(); wu.test.p <- c(); sa.test.p <- c(); moran.st <- c(); moran.p <- c();
for(i in 1:2){
  X <- rayons2[subz.,]
  X <- X[!is.na(as.data.frame(X)[,strsplit(as.character(modz[[i]]$formula),"~")[[2]]]),]
  W_nb <- poly2nb(X, queen = T)
  W_listw <- nb2listw(W_nb, style = "W", zero.policy=T)
  digz.mat["wi.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,3]
  digz.mat["wu.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,3]
  digz.mat["sa.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,3]
  digz.mat["wi.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,4]
  digz.mat["wu.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,4]
  digz.mat["sa.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,4]
  digz.mat["moran.st",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$statistic
  digz.mat["moran.p",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$p.value
}
digz.st <- digz.mat[!grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- digz.mat[grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- ifelse(digz.p<.01,"**",ifelse(digz.p<.05,"*",ifelse(digz.p<.1,"'","")))
for(j in 1:ncol(digz.st)){digz.st[,j] <- as.character(round(digz.st[,j],3))}
for(i in 1:nrow(digz.st)){
  for(j in 1:ncol(digz.st)){
    digz.st[i,j] <- paste0(digz.st[i,j],digz.p[i,j])
  }}
#print.xtable(digz.st,file=gsub("\\.tex","_stats.tex",output.file))
xtable(digz.st)


###########################################################
## Table S10 in Appendix (w/ 2001 census data as covariate)
###########################################################

cv <- "ev_reb_any + ra_agro2 + ra_occupation_duration + ra_urban + ra_ruminukr"
pro.r <- c("po_r14parl_mar","po_r14pres_mar","po_r12parl_mar","po_r10pres_mar","po_r07parl_mar","po_r06parl_mar","po_r04pres_mar")
y.vec <- pro.r
x.vec <- c("ev_deport","ra_ptsn_ops")
z.vec <- c("ra_near_railrd","ra_forest")
X <- rayons2[subz.,]
cv.vec <- strsplit(cv," \\+ ")[[1]]
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
yx.vec <- data.frame(y.vec=rep(y.vec,length(x.vec)),x.vec=rep(x.vec,each=length(y.vec)))
for(j in 1:ncol(yx.vec)){yx.vec[,j] <- as.character(yx.vec[,j])}
yx.vec
summary(X[,cv.vec])

# Model 1
i <- 1
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod1 <- ivreg(formz,data=X);  
summary(mod1,diagnostics=T)

# Model 2
i <- 2
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod2 <- ivreg(formz,data=X);  
summary(mod2,diagnostics=T)


# Model 3
i <- 3
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod3 <- ivreg(formz,data=X);  
summary(mod3,diagnostics=T)

# Model 4
i <- 4
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod4 <- ivreg(formz,data=X);  
summary(mod4,diagnostics=T)

# Model 5
i <- 5
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod5 <- ivreg(formz,data=X);  
summary(mod5,diagnostics=T)

# Model 6
i <- 6
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod6 <- ivreg(formz,data=X);  
summary(mod6,diagnostics=T)

# Model 7
i <- 7
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
for(f in 1:length(cv.vec)){X <- X[!is.na(as.data.frame(X)[,cv.vec[f]]),]}
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod7 <- ivreg(formz,data=X);  
summary(mod7,diagnostics=T)


# Make table (deport)
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
# Diagnostics table (deport)
digz <- c("wi.test","wi.test.p","wu.test","wu.test.p","sa.test","sa.test.p","moran.st","moran.p")
digz.mat <- as.data.frame(matrix(NA,nrow=length(digz),ncol=7))
row.names(digz.mat) <- digz
colnames(digz.mat) <- 1:7
modz <- list(mod1,mod2,mod3,mod4,mod5,mod6,mod7)
wi.test <- c(); wu.test <- c(); sa.test <- c(); wi.test.p <- c(); wu.test.p <- c(); sa.test.p <- c(); moran.st <- c(); moran.p <- c();
for(i in 1:7){
  X <- rayons2[subz.,]
  X <- X[!is.na(as.data.frame(X)[,strsplit(as.character(modz[[i]]$formula),"~")[[2]]]),]
  W_nb <- poly2nb(X, queen = T)
  W_listw <- nb2listw(W_nb, style = "W", zero.policy=T)
  digz.mat["wi.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,3]
  digz.mat["wu.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,3]
  digz.mat["sa.test",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,3]
  digz.mat["wi.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[1,4]
  digz.mat["wu.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[2,4]
  digz.mat["sa.test.p",i] <- summary(modz[[i]],diagnostics=T)$diagnostics[3,4]
  digz.mat["moran.st",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$statistic
  digz.mat["moran.p",i] <- moran.test(x=resid(modz[[i]]),listw=W_listw, zero.policy = T)$p.value
}
digz.st <- digz.mat[!grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- digz.mat[grepl(".p",row.names(digz.mat),fixed=T),]
digz.p <- ifelse(digz.p<.01,"**",ifelse(digz.p<.05,"*",ifelse(digz.p<.1,"'","")))
for(j in 1:ncol(digz.st)){digz.st[,j] <- as.character(round(digz.st[,j],3))}
for(i in 1:nrow(digz.st)){
  for(j in 1:ncol(digz.st)){
    digz.st[i,j] <- paste0(digz.st[i,j],digz.p[i,j])
  }}

xtable(digz.st)



#############################################
## ADDITIONAL IV RESULTS (mentioned in main text, but not reported as tables): RIGHT SECTOR, CPU
#############################################


names(rayons2)

## Regression tables
cv <- "ev_reb_any + ra_agro2 + ra_occupation_duration + ra_urban"
pro.r <- c("po_pravysektor2014", "po_communist12","po_communist07","po_communist06")

y.vec <- pro.r
x.vec <- c("ev_deport","ra_ptsn_ops")
z.vec <- c("ra_near_railrd","ra_forest")
X <- rayons2[subz.,]
cv.vec <- strsplit(cv," \\+ ")[[1]]
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
yx.vec <- data.frame(y.vec=rep(y.vec,length(x.vec)),x.vec=rep(x.vec,each=length(y.vec)))
for(j in 1:ncol(yx.vec)){yx.vec[,j] <- as.character(yx.vec[,j])}
yx.vec

hist(X$ev_reb_any)

# Model 1
i <- 1
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod1 <- ivreg(formz,data=X);  summary(mod1,diagnostics=T)
# first stage
mod1b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod1b,diagnostics=T)


# Model 2
i <- 2
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod2 <- ivreg(formz,data=X);  summary(mod2,diagnostics=T)
# first stage
mod2b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod2b,diagnostics=T)



# Model 3
i <- 3
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod3 <- ivreg(formz,data=X);  
summary(mod3,diagnostics=T)
# first stage
mod3b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod3b,diagnostics=T)

# Model 4
i <- 4
X <- rayons2[subz.,]
X <- X[!is.na(as.data.frame(X)[,yx.vec[i,"y.vec"]]),]
W_nb <- poly2nb(X, queen = T)
W_listw <- nb2listw(W_nb, style = "W", zero.policy = T)
X[[yx.vec[i,"y.vec"]]] <- scale(X[[yx.vec[i,"y.vec"]]]) 
X[[yx.vec[i,"x.vec"]]] <- scale(X[[yx.vec[i,"x.vec"]]])
for(k in 1:length(cv.vec)){X[[cv.vec[k]]] <- scale(X[[cv.vec[k]]])}
for(o in 1:length(z.vec)){X[[z.vec[o]]] <- scale(X[[z.vec[o]]]) }
formz <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast + fitted(ray.ME2)|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)")
formz.base <- paste0(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast|",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast ")
ray.ME2 <- SpatialFiltering(paste(yx.vec[i,"y.vec"],"~",yx.vec[i,"x.vec"],"+",cv,"+ra_oblast ") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
ray.ME1 <- SpatialFiltering(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2)") , data = X, nb = W_nb,style="W", ExactEV = T, zero.policy = T)
mod4 <- ivreg(formz,data=X);  
summary(mod4,diagnostics=T)
# first stage
mod4b <- glm(paste(yx.vec[i,"x.vec"],"~",paste(z.vec,collapse="+"),"+",cv,"+ra_oblast + fitted(ray.ME2) + fitted(ray.ME1)"),data=X); summary(mod4b,diagnostics=T)

stargazer(mod1,mod2,mod3,mod4)
(coefficients(mod2)["ev_deport"]+coefficients(mod3)["ev_deport"]+coefficients(mod4)["ev_deport"])/3



################################################################################
########################Regression Discontinuity Design#########################
################################################################################


## Auxiliary functions used below:
strsplitm <- function(x, w, n) unlist(lapply(strsplit(x, w), function(x) x[n]))

# Prepare data
d <- rayons2
M <- poly2nb(d, queen = FALSE, row.names = d@data$ra_rid)
d$deport <- cut(d$ev_deport, c(-1, 0, 300, 500, 700, Inf), include.lowest = TRUE, labels = 1:5)
id <- rep(1:nrow(d), unlist(lapply(M, length)))
d <- data.frame(id = id, oblast = d$ra_oblast[id], q = d$deport[id], x = d$ev_deport[id], y = d$ev_deport[unlist(M)])


#############################################################
## Figure S2 in Appendix
#############################################################

boxplot(y ~ q, data = d, names = c("0", "(0, 300]", "(300, 500]", "(500, 700]", "> 700"), xlab = "Deportations per district", ylab = "Deportations in neighboring districts", cex.lab = 1.2)

# Partial correlation with neighboring rayons (after adjusting for regional effects), reported in Appendix 5.1
d <- aggregate(cbind(x, y) ~ id + oblast, data = d, FUN = mean)
sqrt(summary(lm(residuals(lm(x ~ oblast, data = d)) ~ y, data = d))$r.squared)

# Pre-treatment balance wrt 1921 census (Appendix 5.2):
# Select precincts within 10 km of the historical borders
data <- rdd_balance
data$diff <- data$rayon_repression - data$neighbor_repression
data$d <- with(data, distance*(diff >= 0) - distance*(diff < 0))/1000
data   <- subset(data, abs(d) < 10)

# Main function to produce balance plots and coefficients for outcome y:
rdout <- function(y, title = ""){
out <- rdrobust(y, data$d, q = 1, p = 0)
out <- paste("Coef. = ", roundr(out$coef[3], 2), " (p = ", roundr(out$pv[3], 2), ")", sep = "")
rdplot(y, data$d, title = paste(title, "\n", out, sep = ""), x.lab = "Distance to rayon border (km)", y.lab = "Outcome")
}


#############################################################
## Figure S3 in Appendix
#############################################################

par(mfrow = c(3, 2))
rdout(y = data$orthdx21/data$habit21, title = "Proportion Orthodox")
rdout(y = data$rc21/data$habit21, title = "Roman Catholics")
rdout(y = data$gc21/data$habit21, title = "Greek Catholics")
rdout(y = data$reljew21/data$habit21, title = "Proportion Jewish")
rdout(y = data$ruthen21/data$habit21, title = "Proportion Ruthenians")
rdout(y = data$polish21/data$habit21, title = "Proportion Polish")


#############################################################
## Figure 3 in main text
#############################################################

# RDD analysis (main text of the paper and Appendix 5.3)
data <- rdd
data$el <- paste(data$year, data$type, sep = "-")
data$diff <- data$rayon_repression - data$neighbor_repression
data$d <- with(data, distance*(diff >= 0) - distance*(diff < 0))/1000

# This function selects pairs historical rayons such that one is at least w standard deviations below the mean and another at least w standard deviations above the mean
s <- function(w) {
  ndata <- subset(data, ((scale(rayon_repression) < -w & scale(neighbor_repression) > w) | (scale(rayon_repression) > w & scale(neighbor_repression) < -w)) & abs(d) < 10)
  ndata$rayon <- factor(strsplitm(as.character(ndata$pair), "\\|", 1))
  ndata
  }

ndata <- s(1)

# Figure 3 of the paper
par(mfrow = c(1, 2))
rdplot(scale(ndata$russian_margin), ndata$d, title = "Distance vs Pro-Russian Vote \n (reduced-form)", x.label = "Distance (in kilometers)", xaxt = 'n', y.lim = c(-1, 1), y.label = "Pro-Russian vote-share (normalized)", col.lines = "black", col.dots="grey50", x.lim = c(-7, 7), cex.lab = 1.2)
axis(1, at = seq(-10, 10, by = 1))
arrows(-3.5, -1, -.1, -1)
text(-3, -.8, "Towards border of \n more repressive district")
arrows(3.5, -1, .1, -1)
text(3, -.8, "Towards border of \n less repressive district")
rdplot(scale(ndata$rayon_repression), ndata$d, title = "Distance vs Deportations  \n (first-stage)", x.label = "Distance (in kilometers)", xaxt = 'n', col.lines = "black", col.dots="grey50", y.label = "District-level deportations (normalized)", x.lim = c(-7, 7), cex.lab = 1.2)
axis(1, at = seq(-10, 10, by = 1))
arrows(-3.5, -1.7, -.1, -1.7)
text(-3, -1.4, "Towards border of \n more repressive district")
arrows(3.5, -1.7, .1, -1.7)
text(3, -1.4, "Towards border of \n less repressive district")
#dev.off()


#############################################################
## Table 3 in main text
#############################################################

# Robust RDD estimators (Lines 1 and 2 of Table 3)
rdrob <- rdrobust(scale(ndata$russian_margin), ndata$d, fuzzy = scale(ndata$rayon_repression))

# The function produces TSLS fuzzy RDD estimates for rayons with w standard deviation contrasts (and clustering of S.E.'s) for lines 3 and 4 of Table 3:
rdfuzzy <- function(w, weights = 0, cluster){
  ndata <- s(w)
  fit <- ivreg(scale(russian_margin) ~ poly(d, 4) + scale(rayon_repression), ~ poly(d, 4) + (d > 0), data = ndata, x = TRUE, weights = weights*abs(diff) + (1-weights)*rep(1, nrow(ndata)))
  coeftest(fit, vcov = cluster.vcov(fit, ndata$rayon, force_posdef=TRUE))[length(coefficients(fit)), ]
}

# Put all results togerther to produce Table 3
x1 <- cbind(rdrob$coef[-1], rdrob$se[-1], rdrob$pv[-1])
x2 <- rdfuzzy(1)[ -3]
x3 <- rdfuzzy(1, weights = 1)[ -3]
ests <- roundr(rbind(x1, x2, x3), 3)
rownames(ests) <- c("1. Bias-correctd ", "2. Bias-corrected with robust errors", "3. TSLS without variance weights", "4. TSLS with variance weights")
print(ests)


#############################################################
## Figure S4 in Appendix
#############################################################

# Function to produce results in Figure S4 of the Appendix:
rdfuzzy <- function(w, cluster){
  ndata <- s(w)
  fit <- ivreg(scale(russian_margin) ~ poly(d, 4) + scale(rayon_repression), ~ poly(d, 4) + (d > 0), data = ndata, x = TRUE, weights = abs(diff))
  x <- coeftest(fit, vcov = cluster.vcov(fit, ndata$rayon, force_posdef=TRUE))[length(coefficients(fit)), ]
  list(b = x[1], ci = c(x[1] - qnorm(0.975)*x[2], x[1] + qnorm(0.975)*x[2]), f = summary(fit, diagnostics = TRUE)$diagnostics[1, 3], n = length(fit$y))
}

w <- seq(.5, 1.5, length = 25)
L <- lapply(w, rdfuzzy)
b <- unlist(lapply(L, function(x) x$b))
ci <- do.call("rbind", lapply(L, function(x) x$ci))
n <- unlist(lapply(L, function(x) x$n))

plot(w, b, pch = 16, xaxt = "n", xlab = "Contrast in repression levels between contiguous districts\n (standard deviations above/below the sample mean).", ylim = c(-.5, 0.13), ylab = "RD estimate")
axis(1, at = seq(.5, 1.5, by = 0.1))
abline(h = 0, lty = 3)
segments(w, ci[,1], w, ci[,2])
text(w,-.45, paste("N", n, sep = " = "), srt = 90, cex = 0.8)


###############################
## Table S6 of Appendix (placebo tests for IV analysis)
###############################

# Prepare data
data <- subset(precinct.sp, ra_oblast%in%c("LVOVSKAYA", "VOLYNSKAYA", "TARNOPOL'SKAYA", "ROVENSKAYA", "DROGOBYCHSKAYA", "STANISLAVSKAYA", "CHERNOVITSKAYA", "ZAKARPATSKAYA"))
data$placebo <- data$ra_oblast=="ZAKARPATSKAYA"
data$el <- paste(data$year, data$type)

# This fits reduced form (precinct-level) regression results for placebo and non-placebo cases that give Table 6 of the Appendix:
b0 <- lmList(r_mar ~ ra_near_railrd | el, data = subset(data, !placebo), na.action = na.omit)
b1 <- lmList(r_mar ~ ra_near_railrd | el, data = subset(data, placebo), na.action = na.omit)
f <- function(out) do.call("rbind", lapply(out, function(x) summary(x)$coefficients[2, c(1, 2, 4)]))
tab <- roundr(cbind(f(b1), f(b0)), 2)
colnames(tab) <- rep(c("Coef.", "S.E.", "p-value"), 2)
rownames(tab) <- rownames(f(b1))
print(tab)



#############################################################
## Figure S5 in Appendix (Caetano exogeneity test)
#############################################################

data <- merged_data

# Figure S5 of Appendix
hist(data$ev_deport, breaks = 50, col = "dark blue", xlab = "Number of deported", main = "")

# Reshape the data to pool across all elections and save to Stata format
sdata <- subset(data, select = c(po_r14parl_mar, po_r14pres_mar, po_r12parl_mar, po_r10pres_mar, po_r07parl_mar, po_r04pres_mar, ev_deport, ev_reb_any, ra_agro2, ra_occupation_duration, ra_oblast, ra_dist_to_ra))
sdata <- reshape(sdata, varying = list(1:6), v.names = "margin", direction = "long")
sdata$margin <- unlist(tapply(sdata$margin, sdata$time, scale))
write.dta(sdata, file = "stata_exogen.dta")

# The remaining analysis is conducted in State using the code in archive file Caetano.zip, which also contains data file in Stata format and its own readme.txt file to explain the contents



##########################################
## Table S8 of Appendix (post-treatment adjustment for population)
##########################################

d <- merged_data
d <- with(merged_data, 
          rbind(
            data.frame(election = 1, oblast = ra_oblast, y = po_r14parl_mar, voters = po_voters14parl, ev_deport = ev_deport),
            data.frame(election = 2, oblast = ra_oblast, y = po_r14pres_mar, voters = po_voters14parl, ev_deport = ev_deport),
            data.frame(election = 3, oblast = ra_oblast, y = po_r12parl_mar, voters = po_voters12, ev_deport = ev_deport),
            data.frame(election = 4, oblast = ra_oblast, y = po_r10pres_mar, voters = po_voters10, ev_deport = ev_deport),
            data.frame(election = 5, oblast = ra_oblast, y = po_r07parl_mar, voters = po_voters07, ev_deport = ev_deport),
            data.frame(election = 6, oblast = ra_oblast, y = po_r06parl_mar, voters = po_voters06, ev_deport = ev_deport),
            data.frame(election = 7, oblast = ra_oblast, y = po_r04pres_mar, voters = po_voters04, ev_deport = ev_deport)
          ))

fit <- NULL
fit[[1]] <- lm(scale(y) ~ factor(election) + factor(oblast) + scale(ev_deport), data = d)
fit[[2]] <- update(fit[[1]], . ~ . + scale(log(voters)))
cluster <- lapply(fit, function(fit) sqrt(diag(cluster.vcov(fit, ~ election))))
stargazer(fit, omit = c("factor", "Constant"), se = cluster)




